{
 "cells": [
  {
   "cell_type": "markdown",
   "source": [
    "In this tutorial we will explore how different combinations of image channels can be defined for image generation in Cue, as well as how to define\n",
    "new SV channels and extend the framework to new sequencing platforms.\n",
    "\n",
    "### Creating new SV signal sets using existing SV channels\n",
    "\n",
    "SV-informative signals are defined in the ```SVSignals``` enum inside ```img/constants```. For example ```SVSignals.RD``` corresponds to the read-depth channel\n",
    "and ```SVSignals.RD_CLIPPED``` is the clipped-read channel.\n",
    "\n",
    "To generate a specific combination of image channels we need to define a new ```SV_SIGNAL_SET``` type and provide its channel composition\n",
    "as a new entry in ```SV_SIGNALS_BY_TYPE``` dictionary inside ```img/constants```. For example, we can create a new set called ```DEMO``` that includes\n",
    "the read-depth and clipped-read signal as follows:\n",
    "\n",
    "```\n",
    "SV_SIGNAL_SET = Enum(\"SV_SIGNAL_SET\", 'SHORT '\n",
    "                                      ...\n",
    "                                      *'DEMO'*)\n",
    "SV_SIGNALS_BY_TYPE = {\n",
    "    SV_SIGNAL_SET.SHORT: [SVSignals.RD, SVSignals.RD_LOW, SVSignals.SR_RP, SVSignals.LLRR, SVSignals.RL,\n",
    "                          SVSignals.LLRR_VS_LR],\n",
    "    ...\n",
    "    *SV_SIGNAL_SET.DEMO: [SVSignals.RD, SVSignals.RD_CLIPPED]*\n",
    "}\n",
    "```\n",
    "Since both of these signal types are already supported by the framework, we just need to specify this new signal set in the data generation YAML\n",
    "config in order to switch images to this channel configuration. In particular, we need to update the ```signal_set``` in ```generate.yaml``` as follows:\n",
    "\n",
    "```\n",
    "signal_set: \"DEMO\"\n",
    "```\n",
    "\n",
    "Providing the updated ```generate.yaml``` config file to the generate.py script will generate 2-channel images with the read-depth and clipped-read channels.\n",
    "\n",
    "### Creating custom SV signals\n",
    "\n",
    "In order to create a new SV channel, we need to extend the ```SVSignals``` definition and add this channel to a signal set (we'll use the ```DEMO``` set).\n",
    "For example, we can define a channel ```RD_MAX``` that will compute the scalar maximum read depth across two loci as follows:\n",
    "\n",
    "```\n",
    "class SVSignals(str, Enum):\n",
    "    RD_MAX = \"RD_MAX\"\n",
    "    ...\n",
    "\n",
    "SV_SIGNALS_BY_TYPE = {\n",
    "    SV_SIGNAL_SET.DEMO: [SVSignals.RD, SVSignals.RD_MAX, SVSignals.RD_CLIPPED]\n",
    "}\n",
    "```\n",
    "\n",
    "We can define the function that this channel should compute inside the ```make_image()``` function of the ```img/datasets``` file by providing\n",
    "a conditional clause for this channel:\n",
    "\n",
    "```\n",
    "if signal == constants.SVSignals.RD_MAX:\n",
    "   counts = self.aln_index.scalar_apply(SVSignals.RD, interval_pair.intervalA, interval_pair.intervalB, op=max)\n",
    "```\n",
    "\n",
    "This definition will apply the ```max``` operator to the read count at each pair of loci represented in the image.\n",
    "\n",
    "### Collecting custom alignment features\n",
    "\n",
    "In order to generate channels from custom alignment features that are not yet supported by the framework, we need to add these features\n",
    "to the BAM index file generated by Cue and used to create the image channels.\n",
    "\n",
    "The ```seq/aln_index.py``` file contains the logic for BAM indexing. In particular, the method ```add_by_signal(...read, signal...)``` implements\n",
    "what/how read alignment properties should be extracted to support a specific signal. The ```read``` input is a pysam ```AlignedSegment``` object\n",
    "which can be queried for various alignment properties of a particular read in the BAM file (this function is called on all the reads in the file).\n",
    "\n",
    "For example, in order to support linked-reads, we construct a channel from the barcodes associated with each read.\n",
    "In particular, the split-molecule channel ```SVSignals.SM``` captures how many barcodes where shared by the reads across\n",
    "each pair of loci in the image. In order to compute the intersection of barcodes across two loci, we store the barcodes observed\n",
    "at each locus as follows:\n",
    "```\n",
    "bin_id = get_bin_id(...read...)\n",
    "if signal == SVSignals.SM :\n",
    "    barcode = read.get_tag('BX')\n",
    "    self.bins[signal][chr_id][bin_id].add(barcode)\n",
    "```\n",
    "\n",
    "In the code above, the read barcode is extracted from the ```BX``` tag in the BAM file and added to the set of barcodes associated with\n",
    "the locus of this read. Once the information is collected in the index, the functions ```scalar_apply()``` or ```intersect()``` provided\n",
    "by the ```AlnIndex``` class (in ```aln_index.py```) can be used implement the desired operation over the values stored in\n",
    "the bins associated with a pair of loci. In the case of the ```SM``` signal, we can use ```intersect()``` to find the intersection of the barcode\n",
    "sets at each pair of loci. Additional features specific to a given sequencing platform can be collected and processed using a similar approach."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}